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There are two usual computational methods for linear (waves and instabilities) problem: eigen- 
value (dispersion relation) solver and initial value solver. In fact, we can introduce an idea of the 
combination of them, i.e., we keep time derivative dt term (and other term if have, e.g., kinetic 
dv term), but transform the linear spatial derivatives dx term to ik, which then can reduce the 
computational dimensions. For example, most (fluid and kinetic) normal mode problems can be 
reduced from treating cumbersome PDEs to treating simple ODEs. Examples for MHD waves, cold 
plasma waves and kinetic Landau damping are given, which show to be extremely simple or even 
may be the simplest method for simulating them. [I don't know whether this idea is new, but it 
seems very interesting and useful. So, I choose making it public] 

I. INTRODUCTION AND BASIC IDEA 

It's well known in plasma physics community that for linear (waves and instabilities) problems we have two usual 
computational methods: eigenvalue (dispersion relation) solver and initial value solver. For the former, we transform 
the linear time derivative d/dt and spatial derivatives V to spectral space using — iu and ik; for the latter, we solve 
the original equations directly. In some simple case, the eigenvalue method can be reduced to analytic tractable form, 
e.g., many well known dispersion relations are this type. However, numerical solutions are always OK (except some 
singularity cases). 

Usually, the eigenvalue method is not intuitive and one needs be good at theoretical derivations; the conventional 
initial value method is complicated in computation and cumbersome in data analysis. Typically, the eigenvalue 
method can give all solutions of the system, while the initial value method can only give the 7 max (most unstable) 
solution. 

Can we combine these two methods? The answer is yes. We can keep time derivative d/dt term (and other term 
if have, e.g., kinetic d/dv term), but transform the linear spatial derivatives <9/<9x term to zk, which is still an 
initial value method, but is solved in half (only spatial not temporal) spectral space, then also has characteristics of 
eigenvalue method. For example, we can highlight the non-7 maa; solutions. 

We will show how to do it with examples for normal mode problems in Sec|n] Since this method has been used 
for linear inhomogeneous eigenvalue and nonlinear problem by previous researchers, we will just give some necessary 
descriptions with citations in Sec |III| 

II. NORMAL MODE PARADIGMS 

In most literatures, normal mode and eigenmode are treated as a same concept since they are very similar. But, 
here, we distinguish them 

Normal mode: homogeneous [2], without boundary conditions; 

Eigen mode: inhomogeneous or with boundary conditions, all possible solutions of the system can be expressed by 
proper sum of eigenmodes. 

Using Half Spectral (HS) method, many normal mode problems can be reduced from PDEs to ODEs. Then can be 
solved extreme easily. While, unfortunately, for eigen mode problems, the equations are in lower-dimensions but still 
PDEs. 
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FIG. 1: RK4 to solve eq.Q, as example for half spectral method, red line is dispersion relation solutions 



A. Simple example 

We construct a simple example to show how to use this method. The original equations 
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gives, cj± = k(u a ± Ub), h = ±h- 

If we assume the initial values fx = yh, the ratio of u± is x and 1 — x respectively, then x — (1 — x) = y =>■ x 
(y + l)/2, then the ratio of the amplitudes for us± is 
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(4) 

which means we can control the amplitudes of each modes of the system exactly by set the proper initial values. 

A 4th order Runge-Kutta simulation of eq. ^ is shown in Figfl] We can see the frequencies and amplitudes of each 
modes (u± = 0.24,0.16, A± = 0.2,0.1) are exact as predict. A small mismatch should be caused by the numerical 
discrete. 



B. ES1D kinetic problem 



Usually, we have two initial value method to simulate kinetic problem (e.g., Landau damping), i.e., Vlasov continuity 
solver and PIC method. 
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FIG. 2: Solve eq.Q for Landau damping 



,X 10 



" k=0.6 



x 10 



14 t=319.98 




FIG. 3: Solve eq.(|5| for beam-plasma instability 



To show that half spectral method is not only for fluid problem, we give a kinetic example. At this subsection, the 
electrostatic ID Landau damping and bump-on-tail simulations are given. 
The linearized equations (ion immobile) are 



which gives 



— = -ikvSf + —E — . 
at m ov 

ikSE — —e I Sfdv. 
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(5) 
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(6) 



contains both normal mode (which is independent with initial value) and ballistic mode (which is brought by initial 
value) or phase mixing (see e.g., [B]). 

The initial distribution function fo can be any form, e.g., Maxwellian gives Landau damping, bump-on-tail gives 
beam-plasma instabilities. 
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Fig{2] shows the simulation of Landau damping. One can find very similar results (especially the fourth panel) from 
Vlasov continuity simulation. However, we should notice an unfavorite recurrence effect P] caused by discrete Av, 
which is also found in half spectral simulation. The recurrence time Tr = 2-njkAv. 

Figj3] is the simulation of bump-on-tail problem. Since wc treat linear problem, the unit of the amplitude can be 
arbitrary large. 

One can find, comparing with continuity solver and PIC method, using half spectral method for Landau damping 
simulation is extreme simple and can be more accurate. Comparing with numerical or analytical dispersion relation 
solver, half spectral method does not need treat troublesome or confusing integral contours. Especially, when the 
initial distribution function is not standard and the dispersion relation is hard to solve, the half spectral simulation 
can gives an reasonable solution for benchmark more complicated codes. 



MHD waves 



For MHD waves, we solve 
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= ik x (<5u x Bq). 
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where B = (0, 0, B ), k = (k sin 9, 0, k cos 6), = jpq/po — ~fkT /m, v 2 A — Bq/ MoPo, v p = w/fc. 

Three solutions are fast mode, slow mode and shear Alfven wave (one can find introductions of them in textbooks, 
e.g., 0]) 
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v A cos 2 6. 



A simulation result is shown in Figj4] Again, we can find the simulation result exactly matches the theoretical 
solutions. The simulation is intuitive. The frequency signal in the second panel is taken from Sp and Su y . If we only 
use 5p signal, the shear Alfven wave solution will vanish. 

If one want to go non-ideal MHD, e.g., including resistivity or using anisotropic pressure 6p\\ ^ <5pj_, but wouldn't 
like to do analytical derivations, then half spectral method is an useful choice: it is very simple, intuitive and can give 
solutions exact enough. 

Since all information for linear perturbation variables is kept in the simulation, we can use them for many more 
deeply analysis, e.g., the polarization and so on. 



D. EM cold plasma waves 



Equations are 



ot m, 



dSE 
dSB 

~dt 



= ic 2 k X SB - <5J/e , 



= ik x SE. 



(9) 



where S3 = J2 S n sos s Sv s - And, Bo = (0, 0, Bo), k = (k sin 9, 0, k cos 9), uj cs — e s Bo/m s and lu ps = n s q 2 /eom s . 

For one ion species, the final dispersion relation can be reduced (with heavy calculations) to a fifth order equation 
for uj 2 (see [10] for details). While, using half spectral simulation, this is very easy. A result is shown in Fig|5] 

If we have more than one ion species or with beams, the dispersion relation can be very headache even though the 
problem seems not that complicated. While, using half spectral method, the problem is indeed still very easy. 
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FIG. 4: Solve eq.Q for ideal MHD waves, red line is dispersion relation solutions 
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FIG. 5: Solve eq.|9| for EM cold plasma waves, red line is dispersion relation solutions 



E. Summary 



For normal mode problem, half spectral method may be the simplest method for simulating them, which is intuitive, 
simple and also exact enough. For analytical difficult problems, this method can not only be an intuitive tool but 
also has practical usages. 
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III. EIGENMODE AND NONLINEAR PROBLEMS 



As claimed at the title, half spectral method can be a general (linear) plasma simulation method. So, we also need 
discuss the cigenmode problem and some nonlinear treatments. While, it is found in the literatures that previous 
researchers have given many examples of this. So, here we just give a short description and mention some citations. 



A. Tearing mode 

Eigenmode problems are similar. We take collisional tearing mode as example here. 

A simulation matches the half spectral idea is given by Lee and Fu[7] (see also citations of that paper). For 
simulation, we solve 



dSB x 

dt 
dSB z 



dt 
d5p 

~dt 



d5p 
~dt = 


dp 
-—bu x - 
dx 


dSu x 
- P dx 


dSu x 


(3 dSp 


B z dSB z 


dt 


2p dx 


p dx 


dSu z 


. 06p 
—ia— 1 





iupbu z , 



= iaB z 5u x 
dB z „ 



1 d 2 5B x a 2 



Rr, 



Rm. dx 2 



-SB Z + ia — -SB X , 
P 



SB X 



(10) 



dx 
dp 



B 



dSu x 



dx 
dSu x 



1 d 2 SB z a 2 



Rr, 



-5u x - 7p 
ox ox 



R m dx 2 
ia^/p5u z . 



bB, 



where parameter are a = kl, /3 = 2/zo??oo/ 'B 2 ^, R m = vaI/v, 7- The normalization unit are B^, p^, va — B 2 ^/ PqPoo, 
Poo, I, t = 1/va- 

The results are very well in that paper[3]. So, we can trust that half spectral method is also good for eigenmode 
problem. 

A bad thing is that, for eigenmode problem, we still need solve PDEs. 



B. Hasegawa-Mima equation 

In fact, an usual way for nonlinear drift wave turbulence by Hasegawa-Mima equation $\ simulation is doing in 
spectral space, which is exact the half spectral idea of this manuscript. A very detailed introduction can be found in 
Waltz's lecture notes [5]. To give a rough impression, one of the equations is shown below 

^kit) = (-iw fc +Jk)4>k(t) + \ ^ S(k- ki - k 2 )Vkk 1 k 2 4>k 1 (t)4>k 2 (t)- (11) 

feifc 2 

We should comment here, for nonlinear problem, we need sum all k modes and the complex conjugate should also 
be kept. 



IV. SUMMARY AND COMMENTS 



In this manuscript, we discussed the idea of half spectral method, and showed that it can be a general simula- 
tion method. However, this idea maybe not new, especially for eigenmode and nonlinear problem, many previous 
researchers used it. I don't know whether this method is new for normal mode problems yet. And the name half 
spectral using here is just for convenient. 

For tokamak (or other problems with strong guide field) simulation, a similar method is called flux tube (e.g., [5]), 
which using the idea of symmetry to reduce dimensions because the physics is mainly along magnetic field line then 
we can take it as one coordinate. Poloidal and toroidal mode numbers m and n are often used in flux tube simulation. 
I haven't checked whether the equations for flux tube simulation are the same for half spectral simulation yet. They 
may be exact the same or at least very similar. 
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But, we can tell that spectral or pseudo spectral method are different from this half spectral method, because 
that they need transform back to real space. For a same simulation, spectral or pseudo spectral method can be 
alternative by other discrete methods (e.g., finite difference) and won't influent the simulation results. While, half 
spectral method is independent on discrete methods. 

This idea is partly inspired by " A 2D Has eg aw a- Mima Model of Electrostatic Drift Wave Turbulence" simulation 
by Deng ZHAO (PKU, 2011) and P239C course project (UCI, 2010) " Cold plasma-warm beam interaction!' given by 
Prof. Liu CHEN. The basic idea is the same. The new here is that we find this method can be generalized. 

If we do not care that whether this idea is new or not, we can find half spectral method is simple, interesting and 
useful. MATLAB codes for this manuscript are given in attached files. 
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